Analyse spatiale et territoriale du logement social
Formation Carthageo-Geoprisme 2021
Introduction
Définir le sujet
Soit le sujet : Logements sociaux et classes sociales dans le Val de Marne en 2017
Définir les “logements sociaux” : Logements HLM ? Logements SRU ?
Définir la notion de “classes sociales” ? : Les CSP ? Un regroupement de CSP ? Les revenus ? la génération
Définir la date : Année 2017 uniquement ? Résultats du RP 2017 (2015-2019) ?
Formuler des questions ou des hypothèses
Qu’elles soient justes ou fausses, les hypothèses permettent de cadrer l’analyse.
CSP et logement social : Les logements sociaux sont réservés aux CSP les plus modestes
Âge et logement social : Les logements sociaux sont réservés aux jeunes ménages actifs
Logement social et territoire : Les logements sociaux sont concentrés dans certains quartiers
Logement social, âge et territoires : Les personnes à la retraite quittent les logements sociaux.
Organiser le travail
Sutout dans le cadre d’un groupe !
Ne collecter que les données utiles pour répondre aux questions posées afin de ne pas être tenté de partir dans toutes les directions
Archiver soigneusement les programmes et les résultats afin de pouvoir reproduire ultérieurement les analyses sur une autre période, un autre territoire.
Ne pas attendre d’avoir accumulé tous les résultats pour les commenter¨ car l’analyse peut suggérer des erreurs ou ouvrir de nouvelles pistes.
Partir des questions et non pas des outils faute de quoi on va trouver des réponses (42 …) sans savoir quelle est la question.
Ne pas confondre les niveaux d’agrégation
Les réponses peuvent varier selon le niveau d’agrégation.
individu / groupe social : “Un ouvrier a plus de chance d’habiter en HLM qu’un cadre”
individu / groupe territorial : “Les cadres et les ouvriers n’habitent pas dans les mêmes quartiers”
groupe social / groupe territorial : “Les cadres évitent les quartiers où il y a beaucoup de logements sociaux”
Données territoriales (RP 2017)
Charger les données statistiques
tab_ind<-readRDS("data2021/94/indiv2017.RDS")
head(tab_ind[,1:5],2)FALSE CANTVILLE NUMMI ACHLR AEMMR AGED
FALSE 1 9401 1 3 9 47
FALSE 2 9401 1 3 9 42
Préparer l’analyse
Soit la relation entre logement en HLM (Y) et CSP du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés
Le choix de la population de référence est important
la sélection ou le regroupement des CSP est également important car il va influer sur les résultats du test.
la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage
Sélectionner les individus et les variables
Puisqu’on travaille sur les logements, on ne garde que les personnes de références des ménages (LPRM=1).
#table(tab_ind$AGEMEN8)
tab_sel<- tab_ind %>%
filter(LPRM == 1) %>%
select(CS1,HLML, IPONDI)
knitr::kable(head(tab_sel,4))| CS1 | HLML | IPONDI |
|---|---|---|
| 3 | 1 | 2.9845995 |
| 4 | 2 | 1.1925306 |
| 8 | 2 | 0.8925770 |
| 2 | 2 | 0.8984924 |
Recoder les modalités
On cherche le code des modalités CS1 ezt HLML dans le fichier des métadonnées
meta<-readRDS("data2021/94/indiv2017_meta.RDS")
metasel <- meta %>% filter(COD_VAR %in% c("CS1", "HLML"))
kable(metasel[,c(1,3,4)])| COD_VAR | COD_MOD | LIB_MOD |
|---|---|---|
| CS1 | 1 | Agriculteurs exploitants |
| CS1 | 2 | Artisans, commerçants et chefs d’entreprise |
| CS1 | 3 | Cadres et professions intellectuelles supérieures |
| CS1 | 4 | Professions Intermédiaires |
| CS1 | 5 | Employés |
| CS1 | 6 | Ouvriers |
| CS1 | 7 | Retraités |
| CS1 | 8 | Autres personnes sans activité professionnelle |
| HLML | 1 | Logement appartenant à un organisme HLM |
| HLML | 2 | Logement n’appartenant pas à un organisme HLM |
| HLML | Z | Hors logement ordinaire |
On recode les modalités des deux variables en regroupant certaines CSP
tab_sel$HLML<-as.factor(tab_sel$HLML)
levels(tab_sel$HLML)<-c("HLM-O","HLM-N")
tab_sel$CS1<-as.factor(tab_sel$CS1)
levels(tab_sel$CS1) <- c("ARCAD","ARCAD","ARCAD","INTER",
"EMPOU","EMPOU","RETRA","INACT")
knitr::kable(head(tab_sel,3))| CS1 | HLML | IPONDI |
|---|---|---|
| ARCAD | HLM-O | 2.984599 |
| INTER | HLM-N | 1.192531 |
| INACT | HLM-N | 0.892577 |
Création du tableau de contingence non pondéré (FAUX)
La solution la plus simple semble être l’instruction table() sad
tab_cont<-table(tab_sel$HLML,tab_sel$CS1)
knitr::kable(addmargins(tab_cont))| ARCAD | INTER | EMPOU | RETRA | INACT | Sum | |
|---|---|---|---|---|---|---|
| HLM-O | 6246 | 11785 | 32802 | 12824 | 4728 | 68385 |
| HLM-N | 46158 | 30804 | 38272 | 34803 | 9865 | 159902 |
| Sum | 52404 | 42589 | 71074 | 47627 | 14593 | 228287 |
Création du tableau de contingence pondéré (JUSTE)
On pondère avec wtd.table() du package questionr.
library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$HLML,tab_sel$CS1,
weights = tab_sel$IPONDI)
knitr::kable(round(addmargins(tab_cont_wtd),0))| ARCAD | INTER | EMPOU | RETRA | INACT | Sum | |
|---|---|---|---|---|---|---|
| HLM-O | 15193 | 28327 | 77954 | 31314 | 10381 | 163169 |
| HLM-N | 123760 | 80991 | 99145 | 98388 | 21225 | 423509 |
| Sum | 138952 | 109318 | 177099 | 129702 | 31607 | 586678 |
Analyser les profils
- Tableau non pondéré … légèrement faux !
| ARCAD | INTER | EMPOU | RETRA | INACT | Ensemble | |
|---|---|---|---|---|---|---|
| HLM-O | 11.9 | 27.7 | 46.2 | 26.9 | 32.4 | 30 |
| HLM-N | 88.1 | 72.3 | 53.8 | 73.1 | 67.6 | 70 |
| Total | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 | 100 |
- Tableau pondéré … juste !
| ARCAD | INTER | EMPOU | RETRA | INACT | Ensemble | |
|---|---|---|---|---|---|---|
| HLM-O | 10.9 | 25.9 | 44 | 24.1 | 32.8 | 27.8 |
| HLM-N | 89.1 | 74.1 | 56 | 75.9 | 67.2 | 72.2 |
| Total | 100.0 | 100.0 | 100 | 100.0 | 100.0 | 100.0 |
Visualiser le tableau de contingence
On choisit l’orientation du tableau et on l’affiche avec plot()
mytable<-wtd.table(tab_sel$CS1,tab_sel$HLML,weights = tab_sel$IPONDI)
plot(mytable, main = "Logements HLM par CSP",
sub = "Source : INSEE - RP 2017",
col=c("lightyellow","lightgreen"))Effectuer un test
Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :
mytest<-chisq.test(mytable)
mytest
Pearson's Chi-squared test
data: mytable
X-squared = 44346, df = 4, p-value < 2.2e-16
Visualisation des résidus
Lorsque la relation est significative, on visualise les cases les plus exceptionnelles avec mosaicplot( …, shade = T)
Conclusion
28% des ménages ordinaires du Val de Marne résident en HLM : Ce chiffre ne tient toutefois pas compte de la population concernée qui peut être plus élevée.
La part du logement en HLM varie bien en fonction de la CSP des actifs : 44% des employés et ouvriers, 26% des intermédiaires, 11% des cadres.
Les inactifs sont également nombreux en HLM : 33% d’entre eux sont dans ce cas, recouvrant des situations diverses.
Une partie des retraités demeurent en HLM : Environ 25% d’entre eux.
Localisation aréale (sf)
Le format sf (spatial features)
La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois
- un tableau de données (l’équivalent du fichier .dbf)
- une géométrie (l’équivalent du fichier .shp)
- une projection (l’équivalent du fichier .prj)
Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou dans d’autres formats standards comme GeoJson, la première tâche consiste donc à les convertir au formt sf afin de pouvoir les utiliser facilement dans R. L’importation se fait à l’aide de l’instruction st_read en indiquant juste le nom du fichier .shp à charger. Les autres fichiers (.dbf ou .proj) seront lus également et intégrés dans l’objet qui hérite de la double classe data.frame et sf.
Etapes de préparation des données
Dans notre exemple, nous allons suivre les étapes suivantes :
- Préparer les données statistiques par IRIS dans un data.frame
- Charger un fonds de carte par IRIS au format sf
- Effectuer une jointure entre les deux fichiers par le code IRIS
- Sauvegarder le résultat
- Agréger les données statistiques et géométriques par commune
- Sauvegarder le résultat.
Préparer les données statistiques
On importe le fichier des individus :
tab_ind<-readRDS("data2021/94/indiv2017.RDS")
head(tab_ind[,1:5],3)FALSE CANTVILLE NUMMI ACHLR AEMMR AGED
FALSE 1 9401 1 3 9 47
FALSE 2 9401 1 3 9 42
FALSE 3 9401 1 3 9 6
Agréger les données
On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :
tab_long<- tab_ind %>%
filter(LPRM == 1)%>%
group_by(IRIS,HLML)%>%
summarise(NB=sum(IPONDI))
knitr::kable(head(tab_long,5),digits=2)| IRIS | HLML | NB |
|---|---|---|
| 940020101 | 1 | 470.39 |
| 940020101 | 2 | 1054.91 |
| 940020102 | 1 | 247.79 |
| 940020102 | 2 | 1651.54 |
| 940020103 | 1 | 419.87 |
Pivoter le tableau
Puis on fait “pivoter” le tableau pour l’obtenir en format large :
tab_large <- tab_long %>% pivot_wider(id_cols = IRIS,
names_from = HLML,
names_prefix = "HLM",
values_from = NB,
values_fill = 0)
knitr::kable(head(tab_large,5),digits=2)| IRIS | HLM1 | HLM2 |
|---|---|---|
| 940020101 | 470.39 | 1054.91 |
| 940020102 | 247.79 | 1651.54 |
| 940020103 | 419.87 | 1030.03 |
| 940020104 | 576.24 | 967.49 |
| 940020105 | 648.40 | 1011.51 |
Ajouter de nouvelles variables
On ajoute de nouvelles variables telles que le nombre total de ménage et le % de ménages en HLM :
tab<- tab_large %>% mutate(TOT = HLM1+HLM2,
HLMpct = 100*HLM1/TOT)
knitr::kable(head(tab,5),digits=2)| IRIS | HLM1 | HLM2 | TOT | HLMpct |
|---|---|---|---|---|
| 940020101 | 470.39 | 1054.91 | 1525.30 | 30.84 |
| 940020102 | 247.79 | 1651.54 | 1899.34 | 13.05 |
| 940020103 | 419.87 | 1030.03 | 1449.90 | 28.96 |
| 940020104 | 576.24 | 967.49 | 1543.73 | 37.33 |
| 940020105 | 648.40 | 1011.51 | 1659.90 | 39.06 |
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par IRIS.
p <- ggplot(tab) + aes (x = HLMpct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90, 100)) +
scale_x_continuous("% de ménages en HLM") +
scale_y_continuous("Nombre d'IRIS") +
ggtitle(label = "Distribution des logements sociaux dans le Val de Marne",
subtitle = "Source : INSEE, RP 2017")
pCharger les données géométriques
On importe le fichier des iris du Val-de-Marne qui est au format sf en ne gardant que les colonnes utiles
map_iris <- readRDS("data2021/94/map_iris.RDS")
map_iris<-map_iris[,c(4,5,1,2,7)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")
class(map_iris)FALSE [1] "sf" "data.frame"
knitr::kable(head(as.data.frame(map_iris)[,1:4],2))| IRIS | NOM_IRIS | COM | NOM_COM |
|---|---|---|---|
| 940210107 | Sorbiers | 94021 | Chevilly-Larue |
| 940680204 | Le Vieux Saint-Maur 4 | 94068 | Saint-Maur-des-Fossés |
Visualisation du fonds iris avec sf
On peut facilement produire une carte vierge des iris du Grand Paris en faisant un plot de la colonne geometry du fichier sf
plot(map_iris$geometry,col="lightyellow")Jointure des données IRIS et du fonds de carte
map_iris_tab<-merge(map_iris,tab,
by.x="IRIS",by.y="IRIS",
all.x=T,all.y=F)
knitr::kable(head(map_iris_tab,3),digits=2)| IRIS | NOM_IRIS | COM | NOM_COM | HLM1 | HLM2 | TOT | HLMpct | geometry |
|---|---|---|---|---|---|---|---|---|
| 940010000 | Ablon-sur-Seine | 94001 | Ablon-sur-Seine | NA | NA | NA | NA | MULTIPOLYGON (((658439.7 68… |
| 940020101 | Chinagora Berthelot | 94002 | Alfortville | 470.39 | 1054.91 | 1525.30 | 30.84 | MULTIPOLYGON (((656616.9 68… |
| 940020102 | Tony Garnier Soladier | 94002 | Alfortville | 247.79 | 1651.54 | 1899.34 | 13.05 | MULTIPOLYGON (((656988.2 68… |
Sauvegarde du fichier par IRIS
On sauvegarde notre fichier au format .RDS de R
saveRDS(map_iris_tab,"data2021/94/map_iris_hlm.RDS")Agrégation statistique + géométriques
Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”
Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.
Agrégation des IRIS en communes
L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries
programme
map_com_tab <- map_iris_tab %>%
group_by(COM, NOM_COM) %>%
summarise(HLM1=sum(HLM1,na.rm=T),
HLM2=sum(HLM2,na.rm=T)) %>%
st_cast("MULTIPOLYGON")
map_com_tab <- map_com_tab %>% mutate(TOT = HLM1+HLM2,
HLMpct = 100*HLM1/TOT)
knitr::kable(head(st_drop_geometry(map_com_tab),3),digits=2)| COM | NOM_COM | HLM1 | HLM2 | TOT | HLMpct |
|---|---|---|---|---|---|
| 94001 | Ablon-sur-Seine | 0.00 | 0.00 | 0.00 | NaN |
| 94002 | Alfortville | 7593.21 | 12562.81 | 20156.02 | 37.67 |
| 94003 | Arcueil | 4505.06 | 5357.96 | 9863.02 | 45.68 |
Agrégation des iris en communes
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par Commune.
programme
p <- ggplot(map_com_tab) + aes (x = HLMpct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90,100)) +
scale_x_continuous("% de ménages en HLM") +
scale_y_continuous("Nombre de communes") +
ggtitle(label = "Distribution des logements sociaux dans le Val de Marne",
subtitle = "Source : INSEE, RP 2017")
pSauvegarde du fichier par commune
On sauvegarde notre fichier au format .RDS de R
saveRDS(map_com_tab,"data2021/94/map_com_hlm.RDS")Visualisation (mapsf)
Le package map_sf
Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.
On trouvera la documentation du package mapsf à l’adresse suivante :
Création d’un template cartographique
Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :
map_iris<-readRDS("data2021/94/map_iris_hlm.RDS")
map_com <-readRDS("data2021/94/map_com_hlm.RDS")tracé d’un fonds de carte vierge
La fonction mf_map() avec le paramètre type = "base"permet de tracer une carte vide
mf_map(map_iris, type = "base")Superposition de couches
On peut toutefois ajouter toute une série de paramètres supplémentaire (col=, border=, lwd=, …) et superposer plusieurs fonds de carte avec le paramètre add = TRUE. L’ajout de la fonction layout permet de rajouter un cadre une légende.
# Trace les Iris avec des paramètres
mf_map(map_iris, type = "base",
col = "lightyellow", border="gray80",lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com, type = "base",
col = NA,border="red",lwd=1,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Val de Marne",
credits = "Sources : IGN et INSEE")Ajout d’un thème
On peut finalement modifier l’ensemble de la carte en lui ajoutant une instruction mf_theme() qui peut reprendre des styles existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”, “candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”, “barcelona”) mais aussi créer ses propres thèmes
#Choix du thème
mf_theme("candy")
mf_map(map_iris, type = "base",
col = "lightyellow", border="gray80",lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="red",lwd=1,
add = TRUE)
mf_layout(title = "Theme candy",
credits = "Sources : IGN et INSEE")Ajout de texte
On peut ajouter une couche de texte avec la fonction mf_label(). Par exemple, on va ajouter à la carte précédente le nom des communes
mf_theme("agolalight")
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray80",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=1,
add = TRUE)
# Ajoute les noms des communes
mf_label(map_com,
var="NOM_COM",
cex=0.4,
col="blue",
overlap = FALSE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et Iris du Val de Marne en 2017",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte de stock
Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.
Dans le package mapsf, on réalise ce type de carte à l’aide de la fonction mf_map()en lui donnant le paramètre type="prop".
On va tenter à titre d’exemple de représenter la distribution du nombre de ménages ordinaires occupant un logement HLM par IRIS :
Carte de stock minimale
Les instructions minimales sont les suivantes :
# Trace les contours des communes
mf_map(x= map_iris,
type = "base")
# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris,
type ="prop",
var = "HLM1",
add=TRUE)Mais le résultat est peu satisfaisant car les cercles sont trop grands. Il faut en pratique toujours effectuer un réglage de ceux-ci avec l’instruction inches=
Carte de stock habillée
mf_theme("agolalight")
mf_map(map_iris, type = "base",
col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="black",lwd=1,add = TRUE)
mf_map(map_iris, var = "HLM1",type = "prop",
inches = 0.1, col = "red",leg_pos = "left",
leg_title = "Nombre de ménages", add=TRUE)
mf_layout(title = "Distribution des logements HLM en 2017",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte choroplèthe
Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.
La fonction du package mapsf adaptée aux variables d’intensité est la fonction mf_map()munie du paramètre type = "choro".
On va prendre l’exemple du nombre de voitures par ménage.
Carte choroplèthe minimale
Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).
# Carte choroplèthe
mf_map(
x = map_iris,
var = "HLMpct",
type = "choro")Carte choroplèthe habillée
On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.
# Choisir les classes et la palette
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens", "Reds"))
# Tracer la carte choroplèthe
mf_map( map_iris, var = "HLMpct",type = "choro",
breaks = mybreaks,pal = mypal,
border="white",col_na = "gray80",
leg_title = "% logements HLM",
leg_val_rnd = 0)
# Ajouter les contours des communes
mf_map(map_com, type = "base", col = NA,
border="black",lwd=1,add = TRUE)
# Ajouter un cadre, un titre et des sources
mf_layout(title = "% de ménages en HLM au RP 2017",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte stock + choroplèthe (1)
On peut combiner les deux modes cartographiques par superposition :
mf_theme("agolalight")
# Choisit les classes
mybreaks = c(0,5,10,20,40,80,100)
# Trace la carte choroplèthe
mf_map(
x = map_iris,
var = "HLMpct",
breaks = mybreaks,
# pal=mypal,
type = "choro",
border="white",
col_na = "gray80",
lwd=0.3,
leg_title = "% ménages",
leg_val_rnd = 0,
)
# Ajoute les cercles proportionnels
mf_map(
x =map_iris,
var = "HLM1",
type = "prop",
inches = 0.06,
col = "red",
leg_pos = "right",
leg_title = "Nb ménages",
add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="black",
lwd=1,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les ménages ordinaires en HLM 2017",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte stock + choroplèthe (2)
Mais on peut aussi utiliser le type prop_choro
mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens", "Reds"))
mf_map(map_iris, type = "base",
col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris, var = c("TOT", "HLMpct"),
inches = 0.08, col_na = "grey", pal=mypal,
breaks = mybreaks, nbreaks = 4, lwd = 0.1,
leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
leg_title = c("nb. ménages", "% HLM"),
add = TRUE)
mf_layout(title = "Les ménages ordinaires en HLM dans le Val de Marne au RP 2017",
frame = TRUE, credits = "Sources : IGN et INSEE")Données spatiales (RPLS, 2020)
Source
Le répertoire des logements locatifs des bailleurs sociaux (RPLS) a pour objectif de dresser l’état global du parc de logements locatifs de ces bailleurs sociaux au 1er janvier d’une année. Il est alimenté par les informations transmises par les bailleurs sociaux. La transmission des informations pour la mise à jour annuelle du répertoire des logements locatifs est obligatoire. Les données sont ensuite géolocalisées à l’adresse et mis à disposition des utilisateurs sur le site du ministère de la transition écologique
Les fichiers sont disponibles en général par régions mais livrés par départements dans le cas de l’Ile de France. Nous allons utilisé ici le fichier du 1er janvier 2020 accessible à l’adresse suivante
https://www.statistiques.developpement-durable.gouv.fr/le-parc-locatif-social-au-1er-janvier-2020-0
Métadonnées
Le fichier de données brutes au format .csv est accompagné d’un document excel précisant le code des variables et la façon dont elles ont été obtenues.
Géocodage
Le fichier indique pour chaque logement sa localisation précise en terme d’adresse mais aussi d’étage dans un immeuble. A partir de ces données qualitatives, l’INSEE a procédé à un géocodage qui aboutit à la création de deux champs :
- coordonnées de latitude et longitude non projetées
- coordonnées de position en projection Lambert officielle
Selon les analyses on peut utiliser l’une ou l’autre de ces coordonnées. Mais la meilleur solution consiste à créer un fichier de type sf (spatial features) en coordonnées WGS94 qu’on pourra ensuite reprojeter dans le système de son choix.
Stratégie
Avant toute exploitation du fichier il est fortement recommandé d’analyser en détail les métadonnées et de définir une stratégie d’analyse.
- choisir une première zone d’étude de petite taille et localisée de préférence dans un espace que l’on connaît bien.
- choisir des variables intéressantes dont l’on connaît bien la signification et dont on a analysé en détail les métadonnées
- vérifier la qualité des données en regardant notamment le nombre de valeurs manquantes, le dégré de précision, etc.
- sélectionner des données auxiliaires issues d’autres sources que l’on souhaite croiser avec celles du RPLS en s’assurant de leur compatibilité (espace, temps, définition, …)
- Ajouter les coordonnées spatiales et stocker le résultat dans un fichier de type sf comportant les indications de projection.
Importation du fichier .csv
Le fichier initial au format .csv est importé à l’aide de la fonction fread() du package data.table car elles est rapide et relativement robuste face aux erreurs de codage. Il faut tout de même préciser le type d’encodage avec encoding ="UTF-8" ainsi que le caractère décimal avec dec=",".
library(data.table)
rpls <- fread("data2021/94/geoloc2020_detail_IDF_dep_94.csv",
encoding = "UTF-8",
dec=",")
rpls<-as.data.frame(rpls)
saveRDS(rpls,"data2021/94/RPLS2020.RDS")Le résultat est converti en data.frame() puis stocké dans un fichier .RDS.
Examen du fichier .RDS
On importe le fichier enregistré au format RDS et on vérifie sa taille avec dim() et sont ype avec class()
don <- readRDS("data2021/94/RPLS2020.RDS")
dim(don)[1] 177902 73
class(don)[1] "data.frame"
Il s’agit bien d’un tableau de données classique de type data.frame. Il comporte 177902 lignes (chacune correspondant à un logement) et 73 variables (décrites dans les métadonnées).
Choix de la zone d’étude
On décide de limiter notre analyse dans un premier temps à quatre communes voisines présentant des profils différents. On peut tester leur profil de respect de la loi SRU en 2019 sur l’application suivante : https://www.ecologie.gouv.fr/sru/
- Bonneuil-sur-Marne (94011): large excédent ( > 25% , pas de pénalité)
- Chennevières-sur-Marne (94019) : léger déficit (22.76%, 58 k€)
- Sucy-en-Brie: déficit (94071) : (19.93% , 150k€ )
- Ormesson-sur Marne (94055) : très fort déficit (2.29%, 665 k€)
On relève leur code INSEE afin de pouvoir faciliter l’extraction des données.
Choix des variables
On va se limiter ici à un très petit nombre de variables
variables de localisation
- result_id : code de l’adresse
- result_label : label de l’adresse
- LIBCOM : nom de la commune
- DEPCOM : code de la commune
- latitude : coordonnées latitude
- longitude: coordonnée longitude
- X : coordonnée projetée (EPSG = 2154)
- Y : coordonnée projetée (EPSG = 2154)
variables thématiques
- CONSTRUCT : année de construction
- SURFHAB : surface habitable en m2
- NBPIECE : nombre de pièces
Extraction du fichier
On applique la double sélection des individus et des variables en nous servant des fonctions filter()et select()du package dplyr.on aboutit ici à un fichier de 8139 lignes et 11 variables.
sel <- don %>%
filter(DEPCOM %in% c("94011","94019",
"94071","94055")) %>%
select(result_id, result_label,
DEPCOM, LIBCOM,
latitude,longitude,X,Y,
CONSTRUCT,SURFHAB,NBPIECE)
dim(sel)[1] 8139 11
Recodage et typage
Certaines variables doivent être recodées ou changées de type afin de faciliter leur exploitation ultérieure par R.
sel$DEPCOM <- as.character(sel$DEPCOM)
sel$LIBCOM <- as.factor(sel$LIBCOM)
sel$PLG_IRIS <- paste(sel$DEPCOM,sel$PLG_IRIS, sep = "")
sel$SURFHAB <- as.numeric(sel$SURFHAB)Résumé rapide
On analyse rapidement les variables thématiques choisies
| CONSTRUCT | SURFHAB | NBPIECE | |
|---|---|---|---|
| Min. :1902 | Min. : 13.00 | Min. :1.000 | |
| 1st Qu.:1966 | 1st Qu.: 55.00 | 1st Qu.:2.000 | |
| Median :1969 | Median : 68.00 | Median :3.000 | |
| Mean :1977 | Mean : 64.58 | Mean :3.176 | |
| 3rd Qu.:1992 | 3rd Qu.: 77.00 | 3rd Qu.:4.000 | |
| Max. :2019 | Max. :172.00 | Max. :6.000 |
Sauvegarde du fichier
On sauvegarde le fichier obtenu au format .RDS afin de garder le formatage des variables :
saveRDS(sel,"data2021/94/sel_logt.RDS")Localisation spatiale
Retour sur sf
Nous revenons sur le package sf (spatial features) que nous avons déjà rencontré au moment de la création de cartes thématiques par IRIS ou communes à l’aide du package mapsf.
Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des logements sociaux au dessus du fonds de carte des IRIS ou communes.
Mais il pourra aussi servir de base à des cartographies dynamiques permettant de placer les points sur des réseaux de rue et plus généralement sur des “tuiles” cartographiques permettant d’effectur des zoom. On utilisera à cet effet d’autres packages comme leaflet ou sa version simplifiée mapview.
Données ponctuelles
Nous reprenons le fichier de localisation établi au chapitre précédent et nous ne conservons que 6 variables:
logt <- readRDS("data2021/94/sel_logt.RDS") %>%
select(adresse=result_id,
X,Y,
date = CONSTRUCT)| adresse | X | Y | date |
|---|---|---|---|
| 94019_0037_00005 | 666779.5 | 6855840 | 1971 |
| 94019_0037_00001 | 666716.7 | 6855829 | 1971 |
| 94019_0037_00001 | 666716.7 | 6855829 | 1971 |
Données IRIS
Nous chargeons par ailleurs le fichier des IRIS en ne gardant que la zone d’étude :
map_iris <- readRDS("data2021/94/map_iris.RDS") %>%
filter(INSEE_COM %in% c("94011","94019",
"94071","94055"))| INSEE_COM | NOM_COM | IRIS | CODE_IRIS | NOM_IRIS | TYP_IRIS | DEPT |
|---|---|---|---|---|---|---|
| 94011 | Bonneuil-sur-Marne | 0103 | 940110103 | Haut Bonneuil | H | 94 |
| 94055 | Ormesson-sur-Marne | 0101 | 940550101 | Nord | H | 94 |
| 94011 | Bonneuil-sur-Marne | 0101 | 940110101 | Zone d’Activite | A | 94 |
Agrégation par commune
Rappel : on peut agréger les géométries d’un fonds sf. Ici on va créer le fonds de carte des communes.
map_com <- map_iris %>% group_by(INSEE_COM,NOM_COM) %>%
summarise() %>%
st_cast("MULTIPOLYGON")Vérification de la projection
Nous savons que les coordonnées X,Y du fichier logement sont projetées en EPS 2154. Mais quelle est la projection de notre fonds IRIS ? S’agit-il de la même ?
st_crs(map_iris)$proj4string[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=44 +lat_2=49 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"
st_crs(2154)$proj4string[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
A priori il s’agit bien de la même (malgré l’inversion de lat_1 et lat_2) de sortte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS
Test de superposition
par(mar=c(0,0,0,0))
#trace les iris
plot(map_iris$geometry,
col="lightyellow", border="gray70",
lwd=0.2)
# trace les communes
plot(map_com$geometry,
col=NA, lwd=1, add=T)
# ajoute les points
points(x=logt$X,
y=logt$Y,
cex=0.2,
col="red",
pch = 16)fichier des adresses
Nous allons maintenant établir un fichier de localisation des adresses en nous servant de l’identifiant unique fourni par l’INSEE.
adr <- logt %>% select(adresse,X,Y) %>%
filter(duplicated(adresse) == F) %>%
filter(is.na(X) ==F,is.na(Y)==F)
dim(adr)[1] 612 3
On constate qu’il n’y a que 652 adresses différentes alors que notre fichier fait état de 8139 logements. Une adresse regroupe donc en moyenne plus de 10 logements (habitat collectif).
Transformation en fichier sf
La transformation de notre fichier initial au format sf est facile à réaliser avec la fonction st_as_sf() du package sf. Mais il faut prendre garde de bien préciser le système de projection si l’on veut pouvoir ensuite l’utiliser.
map_adr <- st_as_sf(adr, coords = c("X","Y"))
st_crs(map_adr)<- 2154
str(map_adr)Classes 'sf' and 'data.frame': 612 obs. of 2 variables:
$ adresse : chr "94019_0037_00005" "94019_0037_00001" "94019_0037_00003" "94019_0037_00007" ...
$ geometry:sfc_POINT of length 612; first list element: 'XY' num 666780 6855840
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
..- attr(*, "names")= chr "adresse"
Agrégation des logements
Notre nouveau fichier sf permet désormais d’effectuer des jointures avec le fichier des logements sociaux. A titre d’exemple on peut désormais compter le nombre de logements par adresse et leur ancienneté moyenne.
logt_by_adr <- logt %>%
group_by(adresse) %>%
summarise(nblog = n(),
datemoy = mean(date))
kable(head(logt_by_adr,10))| adresse | nblog | datemoy |
|---|---|---|
| 31 | 2014 | |
| 94011_0017_00001 | 10 | 1970 |
| 94011_0017_00003 | 10 | 1970 |
| 94011_0017_00005 | 10 | 1970 |
| 94011_0019_00002 | 25 | 1992 |
| 94011_0019_00004 | 24 | 1992 |
| 94011_0022_00001 | 20 | 1966 |
| 94011_0022_00002 | 20 | 1966 |
| 94011_0022_00003 | 20 | 1966 |
| 94011_0022_00004 | 20 | 1966 |
Jointure
On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :
map_logt <- inner_join(logt_by_adr,map_adr) %>% st_as_sf()Cartographie avec mapsf
On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :
mf_theme("agolalight")
mybreaks = c(1900, 1950, 1960, 1970, 1980, 1990, 2000, 2010, 2020)
library(RColorBrewer)
mypal=brewer.pal(n = 8,name = "Spectral")
mf_map(map_iris, type = "base",
col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="black",lwd=1,add = TRUE)
mf_prop_choro( x = map_logt, var = c("nblog", "datemoy"),
inches = 0.08, col_na = "grey", pal=mypal,
breaks = mybreaks, nbreaks = 4, lwd = 0.1,
leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
leg_title = c("nb. logements", "ancienneté"),
add = TRUE)
mf_layout(title = "Les logements sociaux en 2020",
frame = TRUE, credits = "Sources : IGN et RPLS")Sauvegarde des fichiers cartographiques
On sauvegarde nos différents fichiers cartographiques au format sf relatifs à la zone d’étude.
saveRDS(map_com,"data2021/94/sel_map_com.RDS")
saveRDS(map_iris,"data2021/94/sel_map_iris.RDS")
saveRDS(map_logt,"data2021/94/sel_map_logt.RDS")Cartographie dynamique
Statique ou dynamique ?
- Cartographie statique
- production d’images fixes de qualité
- respect strict des règles de la sémiologie graphique
- choix libre d’une projection adaptée (e.g. EPSG 2154)
- production de documents imprimés à finalité normative ou scientifiques
- Cartographie dynamique
- production d’interfaces consultables dans un navigateur.
- modification possible de l’échelle et de l’arrière-plan
- projection imposée par les “tuiles” (EPSG 4326)
- production de documents interactifs à finalité citoyenne ou exploratoire
Packages R de cartographie dynamique
- leaflet : la référence
- Une librairie javascript non liée à un langage (R, Python, html, …)
- Disponible dans R sous forme de package
- Développement constant
- ggmap : l’empire contre attaque
- des outils cartogaphiques utilisant la syntaxe de tidyverse
- impose désormais un lien avec Google
- tmap : une solution hybride
- permet de passer facilement du mode statique au mode dynamique
- mapview : l’équivalent de mapsf
- mis au point par des développeurs allemands
- facilite l’usage de leaflet
- en progrès constant (mais instable)
Préparation des données
On charge les fichiers au format sf et on les transforme en projection WGS94 (EPSG=4326), condition indispensable pour ajouter des “tuiles” dynamiques lors des zoom.
map_com <- readRDS("data2021/94/sel_map_com.RDS") %>%
st_transform(4326)
map_iris <- readRDS("data2021/94/sel_map_iris.RDS") %>%
st_transform(4326)
map_logt <- readRDS("data2021/94/sel_map_logt.RDS") %>%
st_transform(4326)Carte par défaut
Mapview produit par défaut une carte dynamique du fichier sf.
mapview(map_logt)Carte par défaut
On peut zoomer sur la carte, changer les tuiles et faire apparaître des informations sur un point
Superposition de couches
On peut créer des couches et les aditionner avec ‘+’ :
m1 = mapview(map_com, zcol = "NOM_COM")
m2 = mapview(map_logt)
m1+m2Exemple complet
On va essayer de reproduire la carte statique faite avec mapsf
# Carte des communes
map1 <- mapview(map_com, lwd=1, legend= FALSE,
alpha.regions = 0.1)
# Carte des iris
map2 <- mapview(map_iris,lwd = 0.3, label= "NOM_IRIS",
legend= FALSE, alpha.regions = 0)
# Carte des logements
map3 <- mapview(map_logt,
zcol = "datemoy",
at = c(1900,1960, 1970,1980,
1990,2000,2010, 2021),
col.regions = brewer.pal(8, "Spectral"),
cex= "nblog")
map1+map2+map3